
clc
clear
close all

addpath('Data files');
addpath('Regression functions');


Nsim = 1e4; %Paper uses 10000 simulations.

%% Run Main regressions

recLb = 20; %lower bound for amount of recession observations in bootstrap

%Execute script data constructs regression variables
loadRegrVar

%Load macro variables including recession indicators
macroStruc = load('macroVariables.mat');

PMI_RAW = macroStruc.PMI; 
PMI = 1*(PMI_RAW<44.5);
NBER = macroStruc.NBER;  


vec_h = [1 3 6 12];

loopData.mats = matsk;
loopData.q = []; %only needed when using control macro variables
opts.doBoot = 1;
opts.recLb = recLb; 
opts.Nsim = Nsim;
    
for h = vec_h %loop over horizons
            
    hs = strcat('h',num2str(h));
        
    loopData.spreads = regrVar.(hs).spreads;
    loopData.returns = regrVar.(hs).returns;
    loopData.rec = PMI(1:end-h);
    loopData.h = h;
    
    regr.(hs) = loopMatRegr(loopData,opts);
    
end    


%% Main Paper Table 1 

table1= [];

table1(1,:,:)=[regr.h1.mat_bv(:,1)*100 regr.h1.mat_bv(:,2) regr.h1.mat_r2*100 ...
              regr.h1.mat_bvDiff(:,1)*100 regr.h1.mat_bvDiff(:,2) ...
              regr.h1.mat_bvDiff(:,3)*100 regr.h1.mat_bvDiff(:,4) regr.h1.mat_r2Diff*100 ...
              regr.h1.mat_bvSwitch(:,3)*100  regr.h1.mat_bvSwitch(:,4) ];
          
table1(2,:,:)=[regr.h1.mat_tv(:,1) regr.h1.mat_tv(:,2) NaN(numk,1) ...
              regr.h1.mat_tvDiff(:,1) regr.h1.mat_tvDiff(:,2) ...
              regr.h1.mat_tvDiff(:,3) regr.h1.mat_tvDiff(:,4) NaN(numk,1) ...
              regr.h1.mat_tvSwitch(:,3) regr.h1.mat_tvSwitch(:,4)];              

table1(3,:,:)=[regr.h1.mat_tvBoot(:,1) regr.h1.mat_tvBoot(:,2) NaN(numk,1) ...
              regr.h1.mat_tvDiffBoot(:,1) regr.h1.mat_tvDiffBoot(:,2) ...
              regr.h1.mat_tvDiffBoot(:,3) regr.h1.mat_tvDiffBoot(:,4) NaN(numk,1) ...
              regr.h1.mat_tvSwitchBoot(:,3) regr.h1.mat_tvSwitchBoot(:,4)];           
          
fprintf('\n Table 1 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'a ','b ','R2\% ','aExp ','bExp ','aDiff ','bDiff ','R2\% ','aRec ','bRec '); 

for k = (24:12:120);             

    str1 = sprintf(strcat(num2str(k),' %0.2f %0.2f %0.1f %0.2f %0.2f %0.2f %0.2f %0.1f %0.2f %0.2f \n'),table1(1,matsk==k,:)); 
    str2 = sprintf(' (%0.2f) (%0.2f) (%0.3f) (%0.2f) (%0.2f) (%0.2f) (%0.2f) (%0.2f) (%0.2f) (%0.2f) \n',table1(2,matsk==k,:)); 
    str3 = sprintf(' [%0.2f] [%0.2f] [%0.3f] [%0.2f] [%0.2f] [%0.2f] [%0.2f] [%0.2f] [%0.2f] [%0.2f] \n',table1(3,matsk==k,:));
    str4 = sprintf(' \n');
    
    str2 = strrep(str2, '(NaN)', ' ');
    str3 = strrep(str3, '[NaN]', ' ');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end

fprintf('\n \n')

%% Main Paper Table 2:

preWindow = 180;

RiskAvCoef = 5;
  
wLb = 0; 
wUb = 1.5; 


%Discrete yields and returns for asset allocation computations
SPREADS_T = [exp(yields(1:end-1,idxk))-exp(repmat(yields(1:end-1,idx12),1,numk)); NaN(1,numk)];                        
                
logp = -repmat(annMats',T,1).*yields;

rf1 = regrVar.h1.rf;

RETURNS_T = [NaN(1,numk); ...
      exp(logp(1+1:end,idxk_1) - logp(1:end-1,idxk)) - exp(repmat(rf1,1,numk))];   
 
RF_T = [NaN(1,numk); (exp(rf1)-1)*ones(1,numk)];      

PMI_RAW_T = PMI_RAW;

ewmaAlloc = utilGains(PMI_RAW_T,RETURNS_T,SPREADS_T,RF_T,wLb,wUb,'EWMA',preWindow,RiskAvCoef);
tenYrAlloc = utilGains(PMI_RAW_T,RETURNS_T,SPREADS_T,RF_T,wLb,wUb,'10y',preWindow,RiskAvCoef);

tableUtil=[];

tableUtil(1,:,:)=[tenYrAlloc.utilGains(:,1) ...
                  ewmaAlloc.utilGains(:,1) ...
                  tenYrAlloc.utilGains(:,2) ...
                  ewmaAlloc.utilGains(:,2)];
       
tableUtil(2,:,:)=[tenYrAlloc.dmUtilGainsTstat(:,1) ...
                  ewmaAlloc.dmUtilGainsTstat(:,1) ...
                  tenYrAlloc.dmUtilGainsTstat(:,2) ...
                  ewmaAlloc.dmUtilGainsTstat(:,2)];

fprintf('\n Table 2 \n')
fprintf(' %5s %5s %5s %5s \n', ...
                   'Cons Roll','Cons EWMA ','Switch Roll','Swich EWMA'); 

for k = (24:12:120);             

    str1 = sprintf(strcat(num2str(k),' %0.2f   %0.2f   %0.2f   %0.2f \n'),tableUtil(1,matsk==k,:)); 
    str2 = sprintf(' (%0.2f)   (%0.2f)   (%0.2f)   (%0.2f) \n',tableUtil(2,matsk==k,:)); 
    str3 = sprintf(' \n');
    
    str2 = strrep(str2, '(NaN)', ' ');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    
end

fprintf('\n \n')

%% Online Appendix Table 1: VLSTR

c = 44.5;
z = PMI_RAW(1:end-1);

vlstr = mainVlstar(regrVar.h1.returns,regrVar.h1.spreads,z,c);

estGamma = vlstr.gamma;

fprintf('estimated vlstr gamma = %0.1f \n',estGamma)

fX_sixsix  = 1./(1+exp(-exp(6.6)/std(z)*(z-c)));
fX_three  = 1./(1+exp(-exp(3)/std(z)*(z-c)));
fX_two  = 1./(1+exp(-exp(2)/std(z)*(z-c)));
fX_onehalf  = 1./(1+exp(-exp(1.5)/std(z)*(z-c)));

mat_1fX = [1-fX_sixsix 1-fX_three 1-fX_two 1-fX_onehalf];
mat_1fXs = {'sixsix' 'three' 'two' 'onehalf'};

loopData = [];
loopData.spreads = regrVar.h1.spreads;
loopData.returns = regrVar.h1.returns;
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;

opts = [];
opts.doBoot = 0;
opts.recLb = recLb; 
opts.Nsim = Nsim;

for i = 1:size(mat_1fX,2); %loop over horizons
          
    recs = mat_1fXs{i};
        
    loopData.rec = mat_1fX(:,i);
    
    regrFx.(recs) = loopMatRegr(loopData,opts);
       
end    

tableVlstr = [];

tableVlstr(1,:,:)=[regrFx.sixsix.mat_bvDiff(:,2) regrFx.sixsix.mat_bvDiff(:,4) ... 
                 regrFx.three.mat_bvDiff(:,2) regrFx.three.mat_bvDiff(:,4) ...
                 regrFx.two.mat_bvDiff(:,2)  regrFx.two.mat_bvDiff(:,4) ... 
                 regrFx.onehalf.mat_bvDiff(:,2) regrFx.onehalf.mat_bvDiff(:,4) ... 
                 ];
          
      
tableVlstr(2,:,:)=[regrFx.sixsix.mat_tvDiff(:,2) regrFx.sixsix.mat_tvDiff(:,4) ...
                 regrFx.three.mat_tvDiff(:,2) regrFx.three.mat_tvDiff(:,4) ... 
                 regrFx.two.mat_tvDiff(:,2)  regrFx.two.mat_tvDiff(:,4) ... 
                 regrFx.onehalf.mat_tvDiff(:,2)  regrFx.onehalf.mat_tvDiff(:,4)];
               
               
tableVlstr(3,:,:)=[NaN(numk,1) regrFx.sixsix.mat_r2Diff(:,1)*100 ...
                 NaN(numk,1) regrFx.three.mat_r2Diff(:,1)*100 ... 
                 NaN(numk,1)  regrFx.two.mat_r2Diff(:,1)*100 ... 
                 NaN(numk,1)  regrFx.onehalf.mat_r2Diff(:,1)*100];

fprintf('\n Online Table 1: VLSTR \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   '6.6Exp ','6.6Diff ','3Exp ','3Diff ', ... 
                   '2Exp ','2Diff ','1.5Exp ','1.5Diff '); 

for k = (24:12:120);             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  \n'),tableVlstr(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  \n',tableVlstr(2,matsk==k,:)); 
    str3 = sprintf('  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  \n',tableVlstr(3,matsk==k,:));
    str4 = sprintf('  \n');  
 
    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end

fprintf('\n \n')       

%% Main Paper Figure 1: Fitted returns
   
NBER1 = NBER(1:end-1); %lose an obs due to regression
dates1 = dates(1:end-1);

fig1MatIdx = find(matsk == 120);

fig1Names = {'Expected excess returns with switching','Expected excess returns with smooth switching'};

plotHeight = 30;
plotWidth = 30;
subPlotsX = 1;
subPlotsY = 2;   
leftEdge = 1.2;
rightEdge = 0.4;   
topEdge = 2;
bottomEdge = 2;
spaceX = 3;
spaceY = 4;
fontSize = 11; 

subPos = subPlotPos(plotWidth,plotHeight,leftEdge,rightEdge,bottomEdge,topEdge,subPlotsX,subPlotsY,spaceX,spaceY);

%setting the Matlab figure
f = figure('visible','on','name','Data');
clf(f);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [plotWidth plotHeight]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 plotWidth plotHeight]);


for j = 1:2

        ax(j) = axes('position',subPos{j,1},'XGrid','off','XMinorGrid','off','FontSize',fontSize,'Box','on','Layer','top');
        hold on
        
        if j == 1
            p1(1) = plot(dates1,regr.h1.mat_fit(fig1MatIdx,:)'*100);
            p1(2) = plot(dates1,regr.h1.mat_fitSwitch(fig1MatIdx,:)'*100);
            set(p1(1),'Color',[0.5 0.5 0.5],'LineWidth',1.5,'LineStyle',':');
            set(p1(2),'Color',[1 0 0],'LineWidth',0.5)

        elseif j == 2
            p11(1) = plot(dates1,regr.h1.mat_fit(fig1MatIdx,:)'*100);
            p11(2) = plot(dates1,regrFx.two.mat_fitSwitch(fig1MatIdx,:)'*100);
            set(p11(1),'Color',[0.5 0.5 0.5],'LineWidth',1.5,'LineStyle',':');
            set(p11(2),'Color',[1 0 0],'LineWidth',.5)
            
        end
        yMin = min(regr.h1.mat_fitSwitch(fig1MatIdx,:)*100);
        yMax = max(regr.h1.mat_fitSwitch(fig1MatIdx,:)*100);
        
        % For shading above zero
        X  = [dates1',fliplr(dates1)'];
        Y1 = [yMax*NBER1' 0*NBER1'];
        p2 = fill(X,Y1,[0.8 0.8 0.8]);
        set(p2,'edgecolor',[0.8 0.8 0.8]);
        
        % For shading below zero
        Y2 = [yMin*NBER1' 0*NBER1'];
        p3 = fill(X,Y2,[0.8 0.8 0.8]);
        set(p3,'edgecolor',[0.8 0.8 0.8]);

        title(fig1Names{j},'Interpreter','Latex','FontSize',fontSize+4);
        
        hold off
        
        datetick('x','yyyy')
        ylabel('Percent','Interpreter','Latex','FontSize',fontSize+2)
        axis tight
         
        % %force shaded areas to the back of the plot
        uistack(p2, 'bottom')
        uistack(p3, 'bottom')
        
        line([dates1(1) ; dates1(end)],[0 ; 0],'Color','k','LineWidth',0.1,'LineStyle','-');

end

hL = legend(ax(1),[p1(1), p1(2)],{'No switching','With switching'},...
    'Orientation','horizontal','box','off','Interpreter','Latex');
set(hL,'Position',[0.39040316153117 0.96643601040575 0.216425819794804 0.0352564109696279],...
    'Units', 'normalized',...
    'FontSize',14,'Interpreter','Latex');

hLL = legend(ax(2),[p11(1), p11(2)],{'No switching','With smooth switching'},...
    'Orientation','horizontal','box','off','Interpreter','Latex');
set(hLL,'Position', [0.39040316153117 0.47455347511004 0.259365738065072 0.0352564109696283],'Units', 'normalized',...
    'FontSize',14,'Interpreter','Latex');

%% Main Paper Figure 2: Partial R^2 

loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
loopData.returns = regrVar.h1.returns;
loopData.rec = NBER(1:end-1);
loopData.spreads = regrVar.h1.spreads;
regrNBER = loopMatRegr(loopData,opts);

fig2MatIdx = ismember(matsk,[24:12:120]);

figure;
b = bar(annMatsk(fig2MatIdx),100*[regrNBER.mat_r2recDiff(fig2MatIdx) regrNBER.mat_r2expDiff(fig2MatIdx)]);
b(1).FaceColor = 'r';
b(2).FaceColor = 'b';
ylabel('R squared (%)')
xlabel('Maturity in years')
legend({'Recessions','Expansions'},'Location','North')
title('R squared in Recessions and Expansions')    

%% Online Appendix Table 2: Alternative recession indicators

CNAI_RAW=macroStruc.GRO;
CNAI=1*(CNAI_RAW<-0.72);
CNAI(isnan(CNAI_RAW)) = NaN;

ADS_RAW = macroStruc.ADS;
ADS=1*(ADS_RAW<-0.8017);
ADS(isnan(ADS_RAW)) = NaN;

RECPROB_RAW = macroStruc.RECPROB;
RECPROB = 1*(RECPROB_RAW>50);
RECPROB(isnan(RECPROB_RAW)) = NaN;

mat_rec = [PMI ADS RECPROB CNAI NBER];
mat_recs = {'PMI' 'ADS' 'PD' 'CNAI' 'NBER'};

loopData = [];
loopData.spreads = regrVar.h1.spreads;
loopData.returns = regrVar.h1.returns;
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;

opts = [];
opts.doBoot = 1;
opts.recLb = recLb; 
opts.Nsim = Nsim;

for i = 1:size(mat_rec,2); %loop over horizons
           
    recs = mat_recs{i};
        
    loopData.rec = mat_rec(1:end-1,i);
    
    regrRec.(recs) = loopMatRegr(loopData,opts);
       
end    

%Intercepts
tableRecInt = [];

tableRecInt(1,:,:)=[regrRec.PMI.mat_bvDiff(:,1) regrRec.PMI.mat_bvDiff(:,3) ...
                   regrRec.ADS.mat_bvDiff(:,1) regrRec.ADS.mat_bvDiff(:,3) ...
                   regrRec.PD.mat_bvDiff(:,1) regrRec.PD.mat_bvDiff(:,3) ... 
                   regrRec.CNAI.mat_bvDiff(:,1)  regrRec.CNAI.mat_bvDiff(:,3) ... 
                   regrRec.NBER.mat_bvDiff(:,1) regrRec.NBER.mat_bvDiff(:,3)]*100;
          
      
tableRecInt(2,:,:)=[regrRec.PMI.mat_tvDiff(:,1) regrRec.PMI.mat_tvDiff(:,3) ...
                   regrRec.ADS.mat_tvDiff(:,1) regrRec.ADS.mat_tvDiff(:,3) ... 
                   regrRec.PD.mat_tvDiff(:,1) regrRec.PD.mat_tvDiff(:,3) ... 
                   regrRec.CNAI.mat_tvDiff(:,1)  regrRec.CNAI.mat_tvDiff(:,3) ... 
                   regrRec.NBER.mat_tvDiff(:,1) regrRec.NBER.mat_tvDiff(:,3)];
               
               
tableRecInt(3,:,:)=[regrRec.PMI.mat_tvDiffBoot(:,1) regrRec.PMI.mat_tvDiffBoot(:,3) ...
                   regrRec.ADS.mat_tvDiffBoot(:,1) regrRec.ADS.mat_tvDiffBoot(:,3) ... 
                   regrRec.PD.mat_tvDiffBoot(:,1) regrRec.PD.mat_tvDiffBoot(:,3) ... 
                   regrRec.CNAI.mat_tvDiffBoot(:,1)  regrRec.CNAI.mat_tvDiffBoot(:,3) ... 
                   regrRec.NBER.mat_tvDiffBoot(:,1) regrRec.NBER.mat_tvDiffBoot(:,3)];

fprintf('Online Table 2 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'aPMIExp ','aPMIDiff ','ADSExp ','ADSDiff ','PDExp ','PDDiff ', ... 
                   'CNAIExp ','CNAIDiff ','NBERExp ','NBERDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableRecInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableRecInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableRecInt(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end

%Slopes
tableRec = [];

tableRec(1,:,:)=[regrRec.PMI.mat_bvDiff(:,2) regrRec.PMI.mat_bvDiff(:,4) ...
                   regrRec.ADS.mat_bvDiff(:,2) regrRec.ADS.mat_bvDiff(:,4) ...
                   regrRec.PD.mat_bvDiff(:,2) regrRec.PD.mat_bvDiff(:,4) ... 
                   regrRec.CNAI.mat_bvDiff(:,2)  regrRec.CNAI.mat_bvDiff(:,4) ... 
                   regrRec.NBER.mat_bvDiff(:,2) regrRec.NBER.mat_bvDiff(:,4)];
          
      
tableRec(2,:,:)=[regrRec.PMI.mat_tvDiff(:,2) regrRec.PMI.mat_tvDiff(:,4) ...
                   regrRec.ADS.mat_tvDiff(:,2) regrRec.ADS.mat_tvDiff(:,4) ... 
                   regrRec.PD.mat_tvDiff(:,2) regrRec.PD.mat_tvDiff(:,4) ... 
                   regrRec.CNAI.mat_tvDiff(:,2)  regrRec.CNAI.mat_tvDiff(:,4) ... 
                   regrRec.NBER.mat_tvDiff(:,2) regrRec.NBER.mat_tvDiff(:,4)];
               
               
tableRec(3,:,:)=[regrRec.PMI.mat_tvDiffBoot(:,2) regrRec.PMI.mat_tvDiffBoot(:,4) ...
                   regrRec.ADS.mat_tvDiffBoot(:,2) regrRec.ADS.mat_tvDiffBoot(:,4) ... 
                   regrRec.PD.mat_tvDiffBoot(:,2) regrRec.PD.mat_tvDiffBoot(:,4) ... 
                   regrRec.CNAI.mat_tvDiffBoot(:,2)  regrRec.CNAI.mat_tvDiffBoot(:,4) ... 
                   regrRec.NBER.mat_tvDiffBoot(:,2) regrRec.NBER.mat_tvDiffBoot(:,4)];
               
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'bPMIExp ','bPMIDiff ','ADSExp ','ADSDiff ','PDExp ','PDDiff ', ... 
                   'CNAIExp ','CNAIDiff ','NBERExp ','NBERDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableRec(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) \n',tableRec(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableRec(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end

fprintf('\n \n')


%% Online Appendix Table 3: Controlling for Existing Macro Predictors

CF = macroStruc.CF;
CFNAI = macroStruc.GRO;
INF = macroStruc.INF;
LN_update = macroStruc.LN_update;
LN_f1 = macroStruc.LN_f1;

matq = [CFNAI INF LN_update CF];
matqs = {'CFNAI' 'INF' 'LN' 'CPinfl'};

loopData = [];
loopData.spreads = regrVar.h1.spreads;
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI(1:end-1);
loopData.mats = matsk;
loopData.h = 1;

opts = [];
opts.doBoot = 1;
opts.recLb = recLb; 
opts.Nsim = Nsim;

for i = 1:size(matq,2); %loop over horizons
           
    qs = matqs{i};
        
    loopData.q = matq(1:end-1,i);
    
    regrMacro.(qs) = loopMatRegr(loopData,opts);
       
end    

%Intercepts
tableMacroInt = [];

tableMacroInt(1,:,:)=[regr.h1.mat_bvDiff(:,1) regr.h1.mat_bvDiff(:,3) ...
                   regrMacro.CFNAI.mat_bvDiff(:,1) regrMacro.CFNAI.mat_bvDiff(:,3) ...
                   regrMacro.INF.mat_bvDiff(:,1) regrMacro.INF.mat_bvDiff(:,3) ... 
                   regrMacro.LN.mat_bvDiff(:,1)  regrMacro.LN.mat_bvDiff(:,3) ... 
                   regrMacro.CPinfl.mat_bvDiff(:,1) regrMacro.CPinfl.mat_bvDiff(:,3)]*100;
          
      
tableMacroInt(2,:,:)=[regr.h1.mat_tvDiff(:,1) regr.h1.mat_tvDiff(:,3) ...
                   regrMacro.CFNAI.mat_tvDiff(:,1) regrMacro.CFNAI.mat_tvDiff(:,3) ... 
                   regrMacro.INF.mat_tvDiff(:,1) regrMacro.INF.mat_tvDiff(:,3) ... 
                   regrMacro.LN.mat_tvDiff(:,1)  regrMacro.LN.mat_tvDiff(:,3) ... 
                   regrMacro.CPinfl.mat_tvDiff(:,1) regrMacro.CPinfl.mat_tvDiff(:,3)];
               
               
tableMacroInt(3,:,:)=[regr.h1.mat_tvDiffBoot(:,1) regr.h1.mat_tvDiffBoot(:,3) ...
                   regrMacro.CFNAI.mat_tvDiffBoot(:,1) regrMacro.CFNAI.mat_tvDiffBoot(:,3) ... 
                   regrMacro.INF.mat_tvDiffBoot(:,1) regrMacro.INF.mat_tvDiffBoot(:,3) ... 
                   regrMacro.LN.mat_tvDiffBoot(:,1)  regrMacro.LN.mat_tvDiffBoot(:,3) ... 
                   regrMacro.CPinfl.mat_tvDiffBoot(:,1) regrMacro.CPinfl.mat_tvDiffBoot(:,3)];
               
fprintf('\n Online Table 3 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'aExp ','aDiff ','CFNAIExp ','CFNAIDiff ','INFExp ','INFDiff ', ... 
                   'LNExp ','LNDiff ','CPinflExp ','CPinflDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableMacroInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableMacroInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableMacroInt(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end               

%Slopes
tableMacro = [];

tableMacro(1,:,:)=[regr.h1.mat_bvDiff(:,2) regr.h1.mat_bvDiff(:,4) ...
                   regrMacro.CFNAI.mat_bvDiff(:,2) regrMacro.CFNAI.mat_bvDiff(:,4) ...
                   regrMacro.INF.mat_bvDiff(:,2) regrMacro.INF.mat_bvDiff(:,4) ... 
                   regrMacro.LN.mat_bvDiff(:,2)  regrMacro.LN.mat_bvDiff(:,4) ... 
                   regrMacro.CPinfl.mat_bvDiff(:,2) regrMacro.CPinfl.mat_bvDiff(:,4)];
          
      
tableMacro(2,:,:)=[regr.h1.mat_tvDiff(:,2) regr.h1.mat_tvDiff(:,4) ...
                   regrMacro.CFNAI.mat_tvDiff(:,2) regrMacro.CFNAI.mat_tvDiff(:,4) ... 
                   regrMacro.INF.mat_tvDiff(:,2) regrMacro.INF.mat_tvDiff(:,4) ... 
                   regrMacro.LN.mat_tvDiff(:,2)  regrMacro.LN.mat_tvDiff(:,4) ... 
                   regrMacro.CPinfl.mat_tvDiff(:,2) regrMacro.CPinfl.mat_tvDiff(:,4)];
               
               
tableMacro(3,:,:)=[regr.h1.mat_tvDiffBoot(:,2) regr.h1.mat_tvDiffBoot(:,4) ...
                   regrMacro.CFNAI.mat_tvDiffBoot(:,2) regrMacro.CFNAI.mat_tvDiffBoot(:,4) ... 
                   regrMacro.INF.mat_tvDiffBoot(:,2) regrMacro.INF.mat_tvDiffBoot(:,4) ... 
                   regrMacro.LN.mat_tvDiffBoot(:,2)  regrMacro.LN.mat_tvDiffBoot(:,4) ... 
                   regrMacro.CPinfl.mat_tvDiffBoot(:,2) regrMacro.CPinfl.mat_tvDiffBoot(:,4)];
       
               
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'bExp ','bDiff ','CFNAIExp ','CFNAIDiff ','INFExp ','INFDiff ', ... 
                   'LNExp ','LNDiff ','CPinflExp ','CPinflDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableMacro(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableMacro(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableMacro(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end               
                     
fprintf('\n \n')               
               
               
               
               
%% Online Appendix Table 4: Switching Macro Predictors

opts = [];
opts.doBoot = 1;
opts.recLb = recLb; 
opts.Nsim = Nsim;

%CFNAI
loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
CFNAIidxStart = 72;
loopData.returns = regrVar.h1.returns(CFNAIidxStart:end,:); 
loopData.rec = PMI(CFNAIidxStart:end-1);
loopData.spreads = repmat(CFNAI(CFNAIidxStart:end-1)/100,1,numk);
regrSwitchMacro.CFNAI = loopMatRegr(loopData,opts);

%INF 
loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
INFidxStart = 200;
loopData.returns = regrVar.h1.returns(INFidxStart:end,:); 
loopData.rec = PMI(INFidxStart:end-1);
loopData.spreads = repmat(INF(INFidxStart:end-1),1,numk);
regrSwitchMacro.INF = loopMatRegr(loopData,opts);


%LN
loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
loopData.returns = regrVar.h1.returns(1:end,:); 
loopData.rec = PMI(1:end-1);
loopData.spreads = repmat(LN_update(1:end-1),1,numk);
regrSwitchMacro.LN = loopMatRegr(loopData,opts);

%CF
loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
CFidxStart = 126;
loopData.returns = regrVar.h1.returns(CFidxStart:end,:); %missing obs
loopData.rec = PMI(CFidxStart:end-1);
loopData.spreads = repmat(CF(CFidxStart:end-1),1,numk);
regrSwitchMacro.CPinfl = loopMatRegr(loopData,opts);

%LN f1
loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
loopData.returns = regrVar.h1.returns(1:end,:); %missing obs
loopData.rec = PMI(1:end-1);
loopData.spreads = repmat(LN_f1(1:end-1)/100,1,numk);
regrSwitchMacro.LN_f1 = loopMatRegr(loopData,opts);

%Intercepts
tableSwitchMacroInt = [];

tableSwitchMacroInt(1,:,:)=[regrSwitchMacro.CFNAI.mat_bvDiff(:,1) regrSwitchMacro.CFNAI.mat_bvDiff(:,3) ...
                   regrSwitchMacro.INF.mat_bvDiff(:,1) regrSwitchMacro.INF.mat_bvDiff(:,3) ...
                   regrSwitchMacro.LN.mat_bvDiff(:,1)  regrSwitchMacro.LN.mat_bvDiff(:,3) ... 
                   regrSwitchMacro.CPinfl.mat_bvDiff(:,1) regrSwitchMacro.CPinfl.mat_bvDiff(:,3) ...
                   regrSwitchMacro.LN_f1.mat_bvDiff(:,1) regrSwitchMacro.LN_f1.mat_bvDiff(:,3)]*100;
          
tableSwitchMacroInt(2,:,:)=[regrSwitchMacro.CFNAI.mat_tvDiff(:,1) regrSwitchMacro.CFNAI.mat_tvDiff(:,3) ...
                   regrSwitchMacro.INF.mat_tvDiff(:,1) regrSwitchMacro.INF.mat_tvDiff(:,3) ...
                   regrSwitchMacro.LN.mat_tvDiff(:,1)  regrSwitchMacro.LN.mat_tvDiff(:,3) ... 
                   regrSwitchMacro.CPinfl.mat_tvDiff(:,1) regrSwitchMacro.CPinfl.mat_tvDiff(:,3) ...
                   regrSwitchMacro.LN_f1.mat_tvDiff(:,1) regrSwitchMacro.LN_f1.mat_tvDiff(:,3)];              
               
tableSwitchMacroInt(3,:,:)=[regrSwitchMacro.CFNAI.mat_tvDiffBoot(:,1) regrSwitchMacro.CFNAI.mat_tvDiffBoot(:,3) ...
                   regrSwitchMacro.INF.mat_tvDiffBoot(:,1) regrSwitchMacro.INF.mat_tvDiffBoot(:,3) ...
                   regrSwitchMacro.LN.mat_tvDiffBoot(:,1)  regrSwitchMacro.LN.mat_tvDiffBoot(:,3) ... 
                   regrSwitchMacro.CPinfl.mat_tvDiffBoot(:,1) regrSwitchMacro.CPinfl.mat_tvDiffBoot(:,3) ...
                   regrSwitchMacro.LN_f1.mat_tvDiffBoot(:,1) regrSwitchMacro.LN_f1.mat_tvDiffBoot(:,3)];  
               
fprintf('\n Online Table 4 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'CFNAIaExp ','CFNAIaDiff ','INFaExp ','INFaDiff ','LNaExp ','LNaDiff ', ... 
                   'CPinflaExp ','CPinflaDiff ','LNf1aExp ','LNf1aDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableSwitchMacroInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableSwitchMacroInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableSwitchMacroInt(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end               
                
%Slopes
tableSwitchMacro = [];                       

tableSwitchMacro(1,:,:)=[regrSwitchMacro.CFNAI.mat_bvDiff(:,2) regrSwitchMacro.CFNAI.mat_bvDiff(:,4) ...
                   regrSwitchMacro.INF.mat_bvDiff(:,2) regrSwitchMacro.INF.mat_bvDiff(:,4) ...
                   regrSwitchMacro.LN.mat_bvDiff(:,2)  regrSwitchMacro.LN.mat_bvDiff(:,4) ... 
                   regrSwitchMacro.CPinfl.mat_bvDiff(:,2) regrSwitchMacro.CPinfl.mat_bvDiff(:,4) ...
                   regrSwitchMacro.LN_f1.mat_bvDiff(:,2) regrSwitchMacro.LN_f1.mat_bvDiff(:,4)];
          
tableSwitchMacro(2,:,:)=[regrSwitchMacro.CFNAI.mat_tvDiff(:,2) regrSwitchMacro.CFNAI.mat_tvDiff(:,4) ...
                   regrSwitchMacro.INF.mat_tvDiff(:,2) regrSwitchMacro.INF.mat_tvDiff(:,4) ...
                   regrSwitchMacro.LN.mat_tvDiff(:,2)  regrSwitchMacro.LN.mat_tvDiff(:,4) ... 
                   regrSwitchMacro.CPinfl.mat_tvDiff(:,2) regrSwitchMacro.CPinfl.mat_tvDiff(:,4) ...
                   regrSwitchMacro.LN_f1.mat_tvDiff(:,2) regrSwitchMacro.LN_f1.mat_tvDiff(:,4)];              
               
tableSwitchMacro(3,:,:)=[regrSwitchMacro.CFNAI.mat_tvDiffBoot(:,2) regrSwitchMacro.CFNAI.mat_tvDiffBoot(:,4) ...
                   regrSwitchMacro.INF.mat_tvDiffBoot(:,2) regrSwitchMacro.INF.mat_tvDiffBoot(:,4) ...
                   regrSwitchMacro.LN.mat_tvDiffBoot(:,2)  regrSwitchMacro.LN.mat_tvDiffBoot(:,4) ... 
                   regrSwitchMacro.CPinfl.mat_tvDiffBoot(:,2) regrSwitchMacro.CPinfl.mat_tvDiffBoot(:,4) ...
                   regrSwitchMacro.LN_f1.mat_tvDiffBoot(:,2) regrSwitchMacro.LN_f1.mat_tvDiffBoot(:,4)];    
               
tableSwitchMacro(4,:,:)=[NaN(numk,1) regrSwitchMacro.CFNAI.mat_r2Diff(:,1)*100 ...
                        NaN(numk,1) regrSwitchMacro.INF.mat_r2Diff(:,1)*100 ...
                        NaN(numk,1) regrSwitchMacro.LN.mat_r2Diff(:,1)*100 ...
                        NaN(numk,1) regrSwitchMacro.CPinfl.mat_r2Diff(:,1)*100 ...
                        NaN(numk,1) regrSwitchMacro.LN_f1.mat_r2Diff(:,1)*100];

          
               
              
fprintf('  %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s  \n', ...
                   'CFNAIbExp ','CFNAIbDiff ','INFbExp ','INFbDiff ','LNbExp ','LNbDiff ', ... 
                   'CPinflbExp ','CPinflbDiff ','LNf1bExp ','LNf1bDiff  \\'); 

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableSwitchMacro(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) \n',tableSwitchMacro(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableSwitchMacro(3,matsk==k,:));
    str4 = sprintf('  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f \n',tableSwitchMacro(4,matsk==k,:));
    str5 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    fprintf('%s', str5);
    
end


%% Online Appendix Table 5: Alternative Yield-based predictors


opts = [];
opts.doBoot = 1;
opts.recLb = recLb; 
opts.Nsim = Nsim;

%forward spread
forwardSpreads_x12 = 12*regrVar.h1.forwardSpreads;

loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI(1:end-1);
loopData.spreads = forwardSpreads_x12;
regrAltYld.Forward = loopMatRegr(loopData,opts);

%Cochrane-Piazzesi
FBStruc = load('FByields.mat');

FBdates = FBStruc.dates;
FByields = log(1+FBStruc.yields/100);
FBannMatsk =  FBStruc.mat_year';
 
FBCPreturns =  repmat(FBannMatsk(2:end)',T-12,1).*FByields(1:end-12,2:end) ... 
                      -repmat(FBannMatsk(1:end-1)',T-12,1).*FByields(1+12:end,1:end-1) ...
                      -repmat(FByields(1:end-12,1),1,4); 

FBCPlhs = mean(FBCPreturns,2);

FBCPforwards = repmat(FBannMatsk(2:end)',T,1).*FByields(:,2:end) ...
             -repmat(FBannMatsk(1:end-1)',T,1).*FByields(:,1:end-1);

FBCPrhs = [ones(T-12,1) FByields(1:end-12,1) FBCPforwards(1:end-12,:)];

FBgamma = robustOLS(FBCPlhs,FBCPrhs,0,1);

FBCP = [ones(T,1) FByields(:,1) FBCPforwards]*FBgamma;     

loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI(1:end-1);
loopData.spreads = repmat(FBCP(1:end-1),1,numk);
regrAltYld.CP = loopMatRegr(loopData,opts);

%RREL
MA1y = filter(1/12*ones(12,1),1,yields(:,idx12));
RRELlogs = 100*[NaN(12,1); yields(13:end,idx12)-MA1y(12:end-1)]; 

loopData = [];
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;
RRELidxStart = 13;
loopData.returns = regrVar.h1.returns(RRELidxStart:end,:);
loopData.rec = PMI(RRELidxStart:end-1);
loopData.spreads =  repmat(RRELlogs(RRELidxStart:end-1)/100,1,numk);
regrAltYld.RREL = loopMatRegr(loopData,opts);

%Intercepts
tableAltYldInt = [];

tableAltYldInt(1,:,:)=[regrAltYld.Forward.mat_bvDiff(:,1) regrAltYld.Forward.mat_bvDiff(:,3) ...
              regrAltYld.RREL.mat_bvDiff(:,1) regrAltYld.RREL.mat_bvDiff(:,3) ...
              regrAltYld.CP.mat_bvDiff(:,1) regrAltYld.CP.mat_bvDiff(:,3)]*100;
  
tableAltYldInt(2,:,:)=[regrAltYld.Forward.mat_tvDiff(:,1) regrAltYld.Forward.mat_tvDiff(:,3) ...
              regrAltYld.RREL.mat_tvDiff(:,1) regrAltYld.RREL.mat_tvDiff(:,3) ...
              regrAltYld.CP.mat_tvDiff(:,1) regrAltYld.CP.mat_tvDiff(:,3)];          
          
tableAltYldInt(3,:,:)=[regrAltYld.Forward.mat_tvDiffBoot(:,1) regrAltYld.Forward.mat_tvDiffBoot(:,3) ...
              regrAltYld.RREL.mat_tvDiffBoot(:,1) regrAltYld.RREL.mat_tvDiffBoot(:,3) ...
              regrAltYld.CP.mat_tvDiffBoot(:,1) regrAltYld.CP.mat_tvDiffBoot(:,3)];    
          
          
fprintf('\n Online Table 5 \n')          
fprintf(' %5s %5s %5s %5s %5s %5s \n', ...
                   'aExp Fwd ','aDiff Fwd ','aExp RREL ','aDiff RREL ','aExp CP ','aDiff CP '); 

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableAltYldInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) \n',tableAltYldInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  \n',tableAltYldInt(3,matsk==k,:));
    str4 = sprintf('   \n');
    
    str2 = strrep(str2, '(NaN)', ' ');
    str3 = strrep(str3, '[NaN]', ' ');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end

%Slopes
tableAltYld = [];

tableAltYld(1,:,:)=[regrAltYld.Forward.mat_bvDiff(:,2) regrAltYld.Forward.mat_bvDiff(:,4) ...
              regrAltYld.RREL.mat_bvDiff(:,2) regrAltYld.RREL.mat_bvDiff(:,4) ...
              regrAltYld.CP.mat_bvDiff(:,2) regrAltYld.CP.mat_bvDiff(:,4)];
  
tableAltYld(2,:,:)=[regrAltYld.Forward.mat_tvDiff(:,2) regrAltYld.Forward.mat_tvDiff(:,4) ...
              regrAltYld.RREL.mat_tvDiff(:,2) regrAltYld.RREL.mat_tvDiff(:,4) ...
              regrAltYld.CP.mat_tvDiff(:,2) regrAltYld.CP.mat_tvDiff(:,4)];          
          
tableAltYld(3,:,:)=[regrAltYld.Forward.mat_tvDiffBoot(:,2) regrAltYld.Forward.mat_tvDiffBoot(:,4) ...
              regrAltYld.RREL.mat_tvDiffBoot(:,2) regrAltYld.RREL.mat_tvDiffBoot(:,4) ...
              regrAltYld.CP.mat_tvDiffBoot(:,2) regrAltYld.CP.mat_tvDiffBoot(:,4)];    
          
tableAltYld(4,:,:)=[NaN(numk,1) regrAltYld.Forward.mat_r2Diff(:,1)*100 ...
                    NaN(numk,1) regrAltYld.RREL.mat_r2Diff(:,1)*100 ...
                    NaN(numk,1) regrAltYld.CP.mat_r2Diff(:,1)*100];          
          
          
fprintf('  %5s %5s %5s %5s %5s %5s \n', ...
                   'bExp Fwd ','bDiff Fwd ','bExp RREL ','bDiff REL ','bExp CP ','bDiff CP '); 

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableAltYld(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) \n',tableAltYld(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableAltYld(3,matsk==k,:));
    str4 = sprintf('  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f  %0.2f  r2=%0.1f \n',tableAltYld(4,matsk==k,:));
    str5 = sprintf('  \n');
    
    str2 = strrep(str2, '(NaN)', ' ');
    str3 = strrep(str3, '[NaN]', ' ');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    fprintf('%s', str5);
     
end
         
         

%% Online Appendix Table 6: Other short rates for spreads

opts = [];
opts.doBoot = 1;
opts.recLb = recLb; 
opts.Nsim = Nsim;

%PCA
[Wpca,Xpca] = pca(yields(:,mats<=120),'NumComponents',3,'Centered',false);

Xpca2 = -Xpca(1:end-1,2);

loopData = [];
loopData.spreads = repmat(Xpca2,1,numk);
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI(1:end-1);
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;

regrAltSpread.Pca = loopMatRegr(loopData,opts);

%Certificates of deposit
cdStruc = load('cd.mat');

cd = log(1+cdStruc.cd/100);

CdIdxStart = 37;
CdSpreads = yields(CdIdxStart:end-1,idxk)-repmat(cd(CdIdxStart:end-1,1),1,numk); 

loopData = [];
loopData.spreads = CdSpreads; 
loopData.returns = regrVar.h1.returns(CdIdxStart:end,:);
loopData.rec = PMI(CdIdxStart:end-1);
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;

regrAltSpread.Cd = loopMatRegr(loopData,opts);
    
%Eurodollar rates
eurdolStruc = load('eurdol.mat');

eurdol = log(1+eurdolStruc.eurdol/100);

EurdolSpreads = yields(1:end-1,idxk)-repmat(eurdol(1:end-1,1),1,numk);

loopData = [];
loopData.spreads = EurdolSpreads;
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI(1:end-1);
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;

regrAltSpread.Eurdol = loopMatRegr(loopData,opts);

%T-bill from GW
gwStruc = load('gwData.mat');

gwDates = datenum(num2str(gwStruc.yyyymm),'yyyymm');

idxLate = (gwDates >= datenum(1961,6,1)) & (gwDates <= datenum(2016,12,31));

gwTblLate = log(1+gwStruc.tbl(idxLate));

TblSpreads = yields(1:end-1,idxk)-repmat(gwTblLate(1:end-1,1),1,numk);

loopData = [];
loopData.spreads = TblSpreads;
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI(1:end-1);
loopData.mats = matsk;
loopData.q = [];
loopData.h = 1;

regrAltSpread.Tbl = loopMatRegr(loopData,opts);


% Intercepts
tableSpreadInt = [];

tableSpreadInt(1,:,:)=[regr.h1.mat_bvDiff(:,1) regr.h1.mat_bvDiff(:,3) ...
                   regrAltSpread.Pca.mat_bvDiff(:,1) regrAltSpread.Pca.mat_bvDiff(:,3) ...
                   regrAltSpread.Cd.mat_bvDiff(:,1) regrAltSpread.Cd.mat_bvDiff(:,3) ... 
                   regrAltSpread.Eurdol.mat_bvDiff(:,1) regrAltSpread.Eurdol.mat_bvDiff(:,3) ... 
                   regrAltSpread.Tbl.mat_bvDiff(:,1) regrAltSpread.Tbl.mat_bvDiff(:,3)]*100;
          
      
tableSpreadInt(2,:,:)=[regr.h1.mat_tvDiff(:,1) regr.h1.mat_tvDiff(:,3) ...
                   regrAltSpread.Pca.mat_tvDiff(:,1) regrAltSpread.Pca.mat_tvDiff(:,3) ... 
                   regrAltSpread.Cd.mat_tvDiff(:,1) regrAltSpread.Cd.mat_tvDiff(:,3) ... 
                   regrAltSpread.Eurdol.mat_tvDiff(:,1) regrAltSpread.Eurdol.mat_tvDiff(:,3) ... 
                   regrAltSpread.Tbl.mat_tvDiff(:,1) regrAltSpread.Tbl.mat_tvDiff(:,3)];
               
               
tableSpreadInt(3,:,:)=[regr.h1.mat_tvDiffBoot(:,1) regr.h1.mat_tvDiffBoot(:,3) ...
                   regrAltSpread.Pca.mat_tvDiffBoot(:,1) regrAltSpread.Pca.mat_tvDiffBoot(:,3) ... 
                   regrAltSpread.Cd.mat_tvDiffBoot(:,1) regrAltSpread.Cd.mat_tvDiffBoot(:,3) ... 
                   regrAltSpread.Eurdol.mat_tvDiffBoot(:,1) regrAltSpread.Eurdol.mat_tvDiffBoot(:,3) ... 
                   regrAltSpread.Tbl.mat_tvDiffBoot(:,1) regrAltSpread.Tbl.mat_tvDiffBoot(:,3)];
               
fprintf('\n Online Table 6 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'aExp ','aDiff ','PCA2Exp ','PCA2Diff ','CDExp ','CDDiff ', ... 
                   'EurdolExp ','EurdolDiff ','TblExp ','TblDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableSpreadInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableSpreadInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableSpreadInt(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end     

% Slopes
tableSpread = [];

tableSpread(1,:,:)=[regr.h1.mat_bvDiff(:,2) regr.h1.mat_bvDiff(:,4) ...
                   regrAltSpread.Pca.mat_bvDiff(:,2) regrAltSpread.Pca.mat_bvDiff(:,4) ...
                   regrAltSpread.Cd.mat_bvDiff(:,2) regrAltSpread.Cd.mat_bvDiff(:,4) ... 
                   regrAltSpread.Eurdol.mat_bvDiff(:,2) regrAltSpread.Eurdol.mat_bvDiff(:,4) ... 
                   regrAltSpread.Tbl.mat_bvDiff(:,2) regrAltSpread.Tbl.mat_bvDiff(:,4)];
          
      
tableSpread(2,:,:)=[regr.h1.mat_tvDiff(:,2) regr.h1.mat_tvDiff(:,4) ...
                   regrAltSpread.Pca.mat_tvDiff(:,2) regrAltSpread.Pca.mat_tvDiff(:,4) ... 
                   regrAltSpread.Cd.mat_tvDiff(:,2) regrAltSpread.Cd.mat_tvDiff(:,4) ... 
                   regrAltSpread.Eurdol.mat_tvDiff(:,2) regrAltSpread.Eurdol.mat_tvDiff(:,4) ... 
                   regrAltSpread.Tbl.mat_tvDiff(:,2) regrAltSpread.Tbl.mat_tvDiff(:,4)];
               
               
tableSpread(3,:,:)=[regr.h1.mat_tvDiffBoot(:,2) regr.h1.mat_tvDiffBoot(:,4) ...
                   regrAltSpread.Pca.mat_tvDiffBoot(:,2) regrAltSpread.Pca.mat_tvDiffBoot(:,4) ... 
                   regrAltSpread.Cd.mat_tvDiffBoot(:,2) regrAltSpread.Cd.mat_tvDiffBoot(:,4) ... 
                   regrAltSpread.Eurdol.mat_tvDiffBoot(:,2) regrAltSpread.Eurdol.mat_tvDiffBoot(:,4) ... 
                   regrAltSpread.Tbl.mat_tvDiffBoot(:,2) regrAltSpread.Tbl.mat_tvDiffBoot(:,4)];
              
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'bExp ','bDiff ','PCA2Exp ','PCA2Diff ','CDExp ','CDDiff ', ... 
                   'EurdolExp ','EurdolDiff ','TblExp ','TblDiff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableSpread(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableSpread(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableSpread(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end



%% Online Appendix Figure 2: Alternative yield spreads

ofig2MatIdx = find(matsk == 120);

plotSpread = regrVar.h1.spreads(:,ofig2MatIdx);

Xpca2Scaled = ((Xpca2 - mean(Xpca2))/std(Xpca2))*std(plotSpread) + mean(plotSpread);

mat_plotSpread = [Xpca2Scaled [NaN(CdIdxStart-1,1); CdSpreads(:,ofig2MatIdx)] EurdolSpreads(:,ofig2MatIdx) TblSpreads(:,ofig2MatIdx)];

spreadNames = {'2nd PCA','CD','Euro-Dollar','T-Bill'};
spreadLegends = {'2nd PCA','$y_{t,120}-y_{t,3}^{CD}$','$y_{t,120}-y_{t,3}^{ED}$','$y_{t,120}-y_{t,3}^{T-bill}$'};

plotHeight = 25;
plotWidth = 30;
subPlotsX = 2;
subPlotsY = 2;   
leftEdge = 1.2;
rightEdge = 0.4;   
topEdge = 1.5;
bottomEdge = 1.5;
spaceX = 1.5;
spaceY = 1.5;
fontSize = 10.5; 

subPos = subPlotPos(plotWidth,plotHeight,leftEdge,rightEdge,bottomEdge,topEdge,subPlotsX,subPlotsY,spaceX,spaceY);

f = figure('visible','on','name','Data');
clf(f);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [plotWidth plotHeight]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 plotWidth plotHeight]);

factorPlot = 0.0;

j = 0;

newUnits = 'normalized';   

for i = 1:subPlotsY
    for ii = 1:subPlotsX
    
    j = j + 1;   
     
    ax(ii) = axes('position',subPos{i,ii},'XGrid','off','XMinorGrid','off','FontSize',fontSize,'Box','on','Layer','top');
    
    hold on
    
    pj = plot(dates1,[plotSpread mat_plotSpread(:,j)]);

    title(strcat('\fontsize{10.8}',spreadNames{j}))
    
    yMin = min(min([plotSpread mat_plotSpread(:,j)]))*1.1;
    yMax = max(max([plotSpread mat_plotSpread(:,j)]))*1.3;  
   
    % For shading above zero
    X  = [dates1',fliplr(dates1)'];
    Y1 = [yMax*NBER1' 0*NBER1'];
    p2 = fill(X,Y1,[0.8 0.8 0.8]);
    set(p2,'edgecolor',[0.8 0.8 0.8]);  

    % For shading below zero
    Y2 = [yMin*NBER1' 0*NBER1'];
    p3 = fill(X,Y2,[0.8 0.8 0.8]);
    set(p3,'edgecolor',[0.8 0.8 0.8]);

    hold off

    datetick('x','yyyy')

    %force shaded areas to the back of the plot
    uistack(p2, 'bottom')
    uistack(p3, 'bottom')   

    line([dates(1) ; dates(end-1)],[0 ; 0],'Color','k','LineWidth',0.1,'LineStyle','-');

    axis([dates(1) dates(end-1) yMin yMax])   
    
    set(pj(1),'Color','b','LineWidth',1);
    set(pj(2),'Color','k','LineWidth',1);
    
    hL = legend(pj([1 2]),{'$y_{t,120}-y_{t,12}$',spreadLegends{j}},'Orientation','horizontal','box','off','interpreter','latex');
    
    set(hL,'Units', newUnits,'FontSize',12,'Location','North');
    
    end
end


%% Online Appendix Table 7: Recession robustness per decade

datest = datetime(dates,'ConvertFrom','datenum');

matDateIdx = [(year(datest)>=1970 & year(datest)<1980) (year(datest)>=1980 & year(datest)<1990) (year(datest)>=1990 & year(datest)<2000) year(datest)>=2000];

matDateIdxs = {'Seventy' 'Eighty' 'Ninety' 'Thousand'};

PMI1 = PMI(1:end-1);

loopData = [];
loopData.spreads = regrVar.h1.spreads;
loopData.returns = regrVar.h1.returns;
loopData.rec = PMI1;
loopData.mats = matsk;
loopData.h = 1;
loopData.q = [];

for i = 1:size(matDateIdx,2); %loop over horizons
          
    dateIdxs = matDateIdxs{i};
        
    %take out decade
    dateIdx_i = ~(matDateIdx(1:end-1,i).*PMI1(:,1));
       
    opts.doBoot = 1;
    loopData.rec = PMI1(dateIdx_i);
    loopData.spreads = regrVar.h1.spreads(dateIdx_i,:);
    loopData.returns = regrVar.h1.returns(dateIdx_i,:);
    
    regrDecade.(dateIdxs) = loopMatRegr(loopData,opts);
       
end    

%Intercepts
tableDecadeInt = [];
tableDecadeInt(1,:,:)=[regr.h1.mat_bvDiff(:,1) regr.h1.mat_bvDiff(:,3) ...
                   regrDecade.Seventy.mat_bvDiff(:,1) regrDecade.Seventy.mat_bvDiff(:,3) ...
                   regrDecade.Eighty.mat_bvDiff(:,1) regrDecade.Eighty.mat_bvDiff(:,3) ... 
                   regrDecade.Ninety.mat_bvDiff(:,1)  regrDecade.Ninety.mat_bvDiff(:,3) ... 
                   regrDecade.Thousand.mat_bvDiff(:,1) regrDecade.Thousand.mat_bvDiff(:,3)]*100;
          
      
tableDecadeInt(2,:,:)=[regr.h1.mat_tvDiff(:,1) regr.h1.mat_tvDiff(:,3) ...
                   regrDecade.Seventy.mat_tvDiff(:,1) regrDecade.Seventy.mat_tvDiff(:,3) ... 
                   regrDecade.Eighty.mat_tvDiff(:,1) regrDecade.Eighty.mat_tvDiff(:,3) ... 
                   regrDecade.Ninety.mat_tvDiff(:,1)  regrDecade.Ninety.mat_tvDiff(:,3) ... 
                   regrDecade.Thousand.mat_tvDiff(:,1) regrDecade.Thousand.mat_tvDiff(:,3)];
               
               
tableDecadeInt(3,:,:)=[regr.h1.mat_tvDiffBoot(:,1) regr.h1.mat_tvDiffBoot(:,3) ...
                   regrDecade.Seventy.mat_tvDiffBoot(:,1) regrDecade.Seventy.mat_tvDiffBoot(:,3) ... 
                   regrDecade.Eighty.mat_tvDiffBoot(:,1) regrDecade.Eighty.mat_tvDiffBoot(:,3) ... 
                   regrDecade.Ninety.mat_tvDiffBoot(:,1)  regrDecade.Ninety.mat_tvDiffBoot(:,3) ... 
                   regrDecade.Thousand.mat_tvDiffBoot(:,1) regrDecade.Thousand.mat_tvDiffBoot(:,3)];
              
fprintf('\n Online Table 7 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'aExp ','aDiff ','70Exp ','70Diff ','80Exp ','80Diff ', ... 
                   '90Exp ','90Diff ','00Exp ','00Diff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableDecadeInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableDecadeInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableDecadeInt(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end                  
            
% Slopes 
tableDecade = [];

tableDecade(1,:,:)=[regr.h1.mat_bvDiff(:,2) regr.h1.mat_bvDiff(:,4) ...
                   regrDecade.Seventy.mat_bvDiff(:,2) regrDecade.Seventy.mat_bvDiff(:,4) ...
                   regrDecade.Eighty.mat_bvDiff(:,2) regrDecade.Eighty.mat_bvDiff(:,4) ... 
                   regrDecade.Ninety.mat_bvDiff(:,2)  regrDecade.Ninety.mat_bvDiff(:,4) ... 
                   regrDecade.Thousand.mat_bvDiff(:,2) regrDecade.Thousand.mat_bvDiff(:,4)];
          
      
tableDecade(2,:,:)=[regr.h1.mat_tvDiff(:,2) regr.h1.mat_tvDiff(:,4) ...
                   regrDecade.Seventy.mat_tvDiff(:,2) regrDecade.Seventy.mat_tvDiff(:,4) ... 
                   regrDecade.Eighty.mat_tvDiff(:,2) regrDecade.Eighty.mat_tvDiff(:,4) ... 
                   regrDecade.Ninety.mat_tvDiff(:,2)  regrDecade.Ninety.mat_tvDiff(:,4) ... 
                   regrDecade.Thousand.mat_tvDiff(:,2) regrDecade.Thousand.mat_tvDiff(:,4)];
               
               
tableDecade(3,:,:)=[regr.h1.mat_tvDiffBoot(:,2) regr.h1.mat_tvDiffBoot(:,4) ...
                   regrDecade.Seventy.mat_tvDiffBoot(:,2) regrDecade.Seventy.mat_tvDiffBoot(:,4) ... 
                   regrDecade.Eighty.mat_tvDiffBoot(:,2) regrDecade.Eighty.mat_tvDiffBoot(:,4) ... 
                   regrDecade.Ninety.mat_tvDiffBoot(:,2)  regrDecade.Ninety.mat_tvDiffBoot(:,4) ... 
                   regrDecade.Thousand.mat_tvDiffBoot(:,2) regrDecade.Thousand.mat_tvDiffBoot(:,4)];
               
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'bExp ','bDiff ','70Exp ','70Diff ','80Exp ','80Diff ', ... 
                   '90Exp ','90Diff ','00Exp ','00Diff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableDecade(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) (%0.2f)  (%0.2f) \n',tableDecade(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableDecade(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end     

%% Online Appendix Table 8: Long sample

gwDates = datenum(num2str(gwStruc.yyyymm),'yyyymm');

gwTbl = log(1+gwStruc.tbl);
gwLty = log(1+gwStruc.lty);
gwLtr = log(1+gwStruc.ltr);

longNber = gwStruc.nber;

GSStruc = load('GS1.mat'); %Constant mat 1-year from FRED

GS1yield = log(1+GSStruc.GS1/100);

idxNames = {'full','early','late'};
idx.full = (gwDates >= datenum(1926,1,1)) & (gwDates <= datenum(2016,12,31));
idx.early = (gwDates >= datenum(1926,1,1)) & (gwDates <= datenum(1961,5,1));
idx.late = (gwDates >= datenum(1961,6,1)) & (gwDates <= datenum(2016,12,31));

idx.gs1 = (gwDates >= datenum(1953,4,1)) & (gwDates <= datenum(1961,5,1));

% Run regressions
NWlags = 2*1;
blockSize = NWlags;

longTbl = gwTbl;
longTbl(idx.late) = yields(:,mats==12);
longTbl(idx.gs1) = GS1yield;

IbbTable = [];

for i = 1:3
      
   idxi = idx.(idxNames{i}); 
    
   ltri = gwLtr(idxi);
   rfi = 1/12*longTbl(idxi);
   ltyi = gwLty(idxi);
   tbli = longTbl(idxi);
   nberi = longNber(idxi);

   returnsi = ltri(2:end)-rfi(1:end-1);
   spreadi = ltyi(1:end-1)-tbli(1:end-1);
   reci = nberi(1:end-1);
   Ti = length(returnsi);
   Treci = sum(reci);
   
   [bv, sebv, r2] = robustOLS(returnsi,[ones(Ti,1) spreadi],NWlags,1);

   tv = bv./sebv;

   [bvDiff,sebvDiff,r2Diff] = robustOLS(returnsi,[ones(Ti,1) spreadi reci spreadi.*reci],NWlags,1); 

   tvDiff = bvDiff./sebvDiff;

   %Stationary bootstrap
   [bvBoot, sebvBoot] = statBootMinRec( returnsi, [ones(Ti,1) spreadi], ...
                                                                             reci, recLb, Nsim, blockSize );
                                                                                                                                         
   tvBoot = bvBoot./sebvBoot;                                                                     
   
   [bvDiffBoot, sebvDiffBoot] = statBootMinRec( returnsi, [ones(Ti,1) spreadi reci spreadi.*reci], ...
                                                                           reci, recLb, Nsim, blockSize );
                                                                                                                                                  
   tvDiffBoot = bvDiffBoot./sebvDiffBoot;
   
   IbbTable.(idxNames{i}) = [Ti    bv(1)*100  bv(2)     r2*100 bvDiff(1)*100 bvDiff(2)     bvDiff(3)*100 bvDiff(4)     r2Diff*100; ...
                          Treci tv(1)      tv(2)     NaN    tvDiff(1)     tvDiff(2)     tvDiff(3)     tvDiff(4)     NaN; ...
                          NaN   tvBoot(1)  tvBoot(2) NaN    tvDiffBoot(1) tvDiffBoot(2) tvDiffBoot(3) tvDiffBoot(4) NaN];

end

fprintf('\n Online Table 8 \n')     

fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s %5s  \n', ...
                   'T (Trec) ','mu ','theta ','R2\% ','muExp ','thetaExp ','muDiff ','thetaDiff ','R2\% '); 
               
fprintf('1926:1-2016:12  %0.0f  %0.2f  %0.2f  %0.1f  %0.2f  %0.2f  %0.2f  %0.2f  %0.1f  \n',IbbTable.full(1,:)); 
fprintf('                (%0.0f)  (%0.2f) (%0.2f) (%0.3f)   (%0.2f) (%0.2f) (%0.2f) (%0.2f) (%0.3f)  \n',IbbTable.full(2,:));    
fprintf('                [%0.0f]  [%0.2f] [%0.2f] [%0.3f]   [%0.2f] [%0.2f] [%0.2f] [%0.2f] [%0.3f]  \n',IbbTable.full(3,:));    

fprintf('1926:1-1961:05  %0.0f  %0.2f  %0.2f  %0.1f  %0.2f  %0.2f  %0.2f  %0.2f  %0.1f  \n',IbbTable.early(1,:)); 
fprintf('                (%0.0f)  (%0.2f) (%0.2f) (%0.3f)   (%0.2f) (%0.2f) (%0.2f) (%0.2f) (%0.3f)  \n',IbbTable.early(2,:));    
fprintf('                [%0.0f]  [%0.2f] [%0.2f] [%0.3f]   [%0.2f] [%0.2f] [%0.2f] [%0.2f] [%0.3f]  \n',IbbTable.early(3,:));  

fprintf('1961:6-2016:12  %0.0f  %0.2f  %0.2f  %0.1f  %0.2f  %0.2f  %0.2f  %0.2f  %0.1f  \n',IbbTable.late(1,:)); 
fprintf('                (%0.0f)  (%0.2f) (%0.2f) (%0.3f)   (%0.2f) (%0.2f) (%0.2f) (%0.2f) (%0.3f)  \n',IbbTable.late(2,:));    
fprintf('                [%0.0f]  [%0.2f] [%0.2f] [%0.3f]   [%0.2f] [%0.2f] [%0.2f] [%0.2f] [%0.3f]  \n',IbbTable.late(3,:));



%% Online Appendix Table 9: multiple horizons

%Intercepts
tableHorInt = [];

tableHorInt(1,:,:)=[regr.h1.mat_bvDiff(:,1) regr.h1.mat_bvDiff(:,3) ...
                   regr.h3.mat_bvDiff(:,1) regr.h3.mat_bvDiff(:,3) ...
                   regr.h6.mat_bvDiff(:,1) regr.h6.mat_bvDiff(:,3) ... 
                   regr.h12.mat_bvDiff(:,1)  regr.h12.mat_bvDiff(:,3)]*100;
          
      
tableHorInt(2,:,:)=[regr.h1.mat_tvDiff(:,1) regr.h1.mat_tvDiff(:,3) ...
                   regr.h3.mat_tvDiff(:,1) regr.h3.mat_tvDiff(:,3) ... 
                   regr.h6.mat_tvDiff(:,1) regr.h6.mat_tvDiff(:,3) ... 
                   regr.h12.mat_tvDiff(:,1)  regr.h12.mat_tvDiff(:,3)];
               
               
tableHorInt(3,:,:)=[regr.h1.mat_tvDiffBoot(:,1) regr.h1.mat_tvDiffBoot(:,3) ...
                   regr.h3.mat_tvDiffBoot(:,1) regr.h3.mat_tvDiffBoot(:,3) ... 
                   regr.h6.mat_tvDiffBoot(:,1) regr.h6.mat_tvDiffBoot(:,3) ... 
                   regr.h12.mat_tvDiffBoot(:,1)  regr.h12.mat_tvDiffBoot(:,3)];

fprintf('\n Online Table 9 \n')
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'a1Exp ','a1Diff ','3Exp ','3Diff ','6Exp ','6Diff ', ... 
                   '12Exp ','12Diff ');

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableHorInt(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) \n',tableHorInt(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableHorInt(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end


%Slopes
tableHor = [];

tableHor(1,:,:)=[regr.h1.mat_bvDiff(:,2) regr.h1.mat_bvDiff(:,4) ...
                   regr.h3.mat_bvDiff(:,2) regr.h3.mat_bvDiff(:,4) ...
                   regr.h6.mat_bvDiff(:,2) regr.h6.mat_bvDiff(:,4) ... 
                   regr.h12.mat_bvDiff(:,2)  regr.h12.mat_bvDiff(:,4)];
          
      
tableHor(2,:,:)=[regr.h1.mat_tvDiff(:,2) regr.h1.mat_tvDiff(:,4) ...
                   regr.h3.mat_tvDiff(:,2) regr.h3.mat_tvDiff(:,4) ... 
                   regr.h6.mat_tvDiff(:,2) regr.h6.mat_tvDiff(:,4) ... 
                   regr.h12.mat_tvDiff(:,2)  regr.h12.mat_tvDiff(:,4)];
               
               
tableHor(3,:,:)=[regr.h1.mat_tvDiffBoot(:,2) regr.h1.mat_tvDiffBoot(:,4) ...
                   regr.h3.mat_tvDiffBoot(:,2) regr.h3.mat_tvDiffBoot(:,4) ... 
                   regr.h6.mat_tvDiffBoot(:,2) regr.h6.mat_tvDiffBoot(:,4) ... 
                   regr.h12.mat_tvDiffBoot(:,2)  regr.h12.mat_tvDiffBoot(:,4)];
               
fprintf(' %5s %5s %5s %5s %5s %5s %5s %5s \n', ...
                   'b1Exp ','b1Diff ','3Exp ','3Diff ','bExp ','6Diff ', ... 
                   '12Exp ','12Diff '); 

for k = [24 60 84 120];             

    str1 = sprintf(strcat(num2str(k),'  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f  %0.2f \n'),tableHor(1,matsk==k,:)); 
    str2 = sprintf('  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f)  (%0.2f) \n',tableHor(2,matsk==k,:)); 
    str3 = sprintf('  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f]  [%0.2f] \n',tableHor(3,matsk==k,:));
    str4 = sprintf('  \n');

    fprintf('%s', str1);
    fprintf('%s', str2);
    fprintf('%s', str3);
    fprintf('%s', str4);
    
end

fprintf('\n \n')


%% Online Appendix Figure 3: yields per decade

FFRstruc = load('FFR.mat');

FFR = FFRstruc.FFR/100;

matNberDateIdx = [dates>=datenum('30-Jan-1970') & dates<=datenum('30-Nov-1970') ...
                   dates>=datenum('31-Dec-1973') & dates<=datenum('31-Mar-1975') ...
                   dates>=datenum('29-Feb-1980') & dates<=datenum('31-Jul-1980') ...
                   dates>=datenum('31-Aug-1981') & dates<=datenum('30-Nov-1982') ...
                   dates>=datenum('31-Aug-1990') & dates<=datenum('28-Mar-1991') ...
                   dates>=datenum('30-Apr-2001') & dates<=datenum('30-Nov-2001') ...
                   dates>=datenum('31-Jan-2008') & dates<=datenum('30-Jun-2009')];
           
% mat_dateIdxs = {'Seventy1' 'Seventy2' 'Eighty1' 'Eighty2' 'Ninety' 'Thousand1' 'Thousand2'};
matNberDateStr = {'1970' '1973' '1980' '1981' '1990' '2001' '2008'};

plotHeight = 25;
plotWidth = 30;
subPlotsX=3;
subPlotsY=3;   
leftEdge=1.2;
rightEdge=0.4;   
topEdge=1.5;
bottomEdge=1.5;
spaceX=1.5;
spaceY=1.5+0.5;
fontSize=10.5;

subPos=subPlotPos(plotWidth,plotHeight,leftEdge,rightEdge,bottomEdge,topEdge,subPlotsX,subPlotsY,spaceX,spaceY);

%setting the Matlab figure
f=figure('visible','on','name','Understand Excess Returns');
clf(f);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [plotWidth plotHeight]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 plotWidth plotHeight]);

j = 0;

%loop to create axes
for i=1:subPlotsX
    for ii=1:subPlotsY
    
    j = j + 1;
    
    if j<=size(matNberDateIdx,2)
    dateIdxIso_j = (matNberDateIdx(1:end,j))==1;
    
    dates_i = dates(dateIdxIso_j);
    
    plotYields_j = [FFR(dateIdxIso_j) yields(dateIdxIso_j,annMats==2) yields(dateIdxIso_j,annMats==5) yields(dateIdxIso_j,annMats==10)];
    
    ax(i)=axes('position',subPos{i,ii},'XGrid','off','XMinorGrid','off','FontSize',fontSize,'Box','on','Layer','top');
    
    p_i = plot(dates_i,plotYields_j*100);
    datetick('x','mm')
    thdl_i = title(matNberDateStr{j},'FontSize',10.8);
    
    set(p_i(1),'LineWidth',2.25); 
   
    axis tight
    
    end
    
    end

end

hL = legend(p_i,{'FFR','2y','5y','10y'},'Orientation','horizontal','box','off');

set(hL,'Position', [0.3 0.97 0.42 0.04],'Units', 'normalized','FontSize',10); 

%% Online Appendix Table 11: The Aggressiveness of the Fed across Recessions and Expansions
load DTARG
rec=NBER(56:end,1);
for t=1:size(rec,1)
if rec(t)==1 
   DTARGrec(:,t)=DTARG(t);   
end
end
for t=1:size(rec,1)
if rec(t)==0 
   DTARGexp(:,t)=DTARG(t);   
end
end

fprintf('\n Online Table 11 \n')
fprintf('\n Average change in recessions \n')
fprintf(' Tightening    Easing\n')
disp([ mean(DTARGrec( DTARGrec>0 )) mean(DTARGrec( DTARGrec<0 )) ])
fprintf('\n Average change in expansions \n')
fprintf(' Tightening    Easing\n')
disp([ mean(DTARGexp( DTARGexp>0 )) mean(DTARGexp( DTARGexp<0 )) ])






